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SUMMARY 

During  the  past  several  years.  Systems,  Science  and 
Software  (S3)  personnel  have  been  actively  engaged  in  a com- 
prehensive program  involving  computer  modeling  of  the  non- 
linear processes  that  characterize  the  near  source  ground 
motion  from  underground  nuclear  explosions,  propagation  of 
the  resultant  stress  waves  through  the  crust  ^nd  upper 
mantle  and  prediction  of  the  ground  motion  recorded  at  tele- 
seismic  distances  from  the  explosion  [Cherry,  et  al.  (1973), 
Bache,  et  al.  (1975)]. 

The  objective  of  the  research  presented  in  this  report 
is  to  determine  the  dependence  of  teleseismic  magnitudes  on 
the  nonlinear  behavior  of  the  near  source  rock  environment. 
Such  information  enables  one  to  express  uncertainties  in  ex- 
plosive yield  estimates  in  terms  of  uncertainties  in  the 
material  properties  and  provides  insight  concerning  the  re- 
quirements for  collection  of  geophysical  data  at  a specific 
test  site. 

A series  of  spherically  symmetric  calculations  have 
been  performed  in  which  the  properties  of  the  near  source 
material  were  varied  systematically.  Output  from  these  cal- 
culations enables  one  to  determine  the  relative  teleseismic 
coupling  efficiency  of  a given  near  source  material.  Also, 
since  seismic  magnitudes  are  directly  related  to  explosive 
yield,  this  parameter  study  permits  an  assessment  of  uncer- 
tainties in  the  estimated  yield  in  terms  of  uncertainties  of 
the  material  properties  of  the  test  site. 
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I . INTRODUCTION 

1.1  OBJECTIVES 

This  report  describes  the  results  of  a theoretical 
study  directed  at  the  determination  of  the  sensitivity  of  the 
equivalent  elastic  source  from  a nuclear  explosion  to  the 
nonlinear  behavior  of  the  near  source  rock  environment.  Un- 
certainties in  y<  eld  estimates  based  on  measured  seismic 
magnitudes  may  te  expressed  in  terms  of  uncertainties  in  the 
near  source  properties.  The  material  parameters  and  the 
aspects  of  the  constitutive  model  which  have  the  greatest 
effect  on  the  equivalent  source  are  identified  and  guidelines 
are  established  for  the  collection  of  geophysical  data  in  the 
near  source  region. 


1.2  PROCEDURE 

Equivalent  elastic  source  calculations  place  severe 
requirements  on  both  the  constitutive  model  and  the  stress 
wave  simulation  code  which  contains  the  model.  The  constitu- 
tive model  must  account  for  the  various  modes  of  material  be- 
havior which  occur  between  the  cavity  boundary  and  the  dis- 
tance at  which  the  material  response  is  elastic.  The  stress 
wave  code  must  not  only  be  able  to  easily  accept  the  constitu- 
tive model  but  also  accurately  propagate  the  stress  wave  into 
the  elastic  region.  The  equivalent  elastic  source  is  then 
obtained  from  the  properties  of  the  displacement  field  in  the 
elastic  region. 

Constitutive  models,  especially  that  portion  of  the 
model  which  pertains  to  the  material  strength,  can  be  fairly 
complex.  In  order  to  isolate  the  effect  of  the  constitutive 
model  on  seismic  coupling  we  first  assumed  that  only  plastic 
flow  would  result  when  the  allowable  material  strength  was 
exceeded  and  an  extensive  parameter  study  was  completed  for 
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a water-tul f rock  environment  containing  various  amounts  of 
air -filled  porosity.  We  then  modified  the  strength  portion 
of  the  constitutive  model  to  simulate  tension  failure,  shear 
failure  and  the  role  of  pore  fluid  pressure  on  effective 
stress. 

Figure  1.1  depicts,  ii  flow  chart  form,  the  organiza- 
tion of  the  parameter  study.  The  material  selected  consists 
of  a mixture  of  tuff  and  water  which  may  contain  air-filled 
voids.  The  properties  of  the  water  and  grain  density  tuff 
constituents  were  held  constant  for  all  calculations.  The 
mass  fraction  of  water  in  the  mixture,  f^,  is  a natural  para- 
meter to  select  for  initial  variation  and  this  choice  is  in- 
dicated in  Fig.  1.1;  the  air-filled  volume  fraction,  4)  , was 

o 

set  to  zero  for  this  set  of  calculations.  The  equation  of 
state  for  the  tuff-water  mixture  is  described  in  Section  2.1 
of  this  report. 

The  remaining  material  parameters  are  varied  one  at  a 

time  while  holding  the  water  mass  fraction  constant  at  0.17. 

Four  of  the  material  parameters  are  directly  related  to  the 

porosity  model  which  is  discussed  in  Section  2.1:  (1)  the 

initial  air-filled  volume  fraction,  <$  ; (2)  the  elastic 

o 

limit  of  the  pressure,  P ; (3)  the  crush  or  compaction  pres- 
sure, P ; (4)  the  bulk  modulus,  k . The  shear  modulus,  y, 
c o 

is  selected  as  one  of  the  independent  elastic  properties  of 

the  porous  mixture.  Given  k , y and  the  ambient  density,  p , 

o o 

the  elastic  P wave  velocity,  c,  is  given  by 


while  the  shear  wave  velocity,  6,  is 


n Pressure 
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Additional  parameters  are  contained  in  the  material 

strength  model  which  is  discussed  in  Section  2.2.  As  noted 

earlier  and  shown  in  Fig.  1.1,  only  plastic  flow  was  permitted 

during  these  calculations.  The  strength  parameters  included 

in  this  sensitivity  study  are:  (1)  the  zero  pressure  yield 

strength,  Y ; (2)  the  maximum  yield  strength,  Y ; (3)  the 
o m 

pressure  at  which  the  yield  strength  reaches  its  maximum 

value,  P . 

m 

One  last  parameter,  the  overburden  pressure,  P , is 

o 

included  in  the  one-dimensional  sensitivity  study  although 
it  depends  on  the  depth  of  burial  and  is  not  a material 
property. 

A detailed  description  of  the  material  models  employed 
in  our  near  source  calculations  is  presented  in  the  next  sec- 
tion. The  results  of  the  sensitivity  study  are  discussed  in 
Section  III  and  conclusions  based  on  those  results  are  sum- 
marized in  Section  IV.  In  Section  V the  strength  portion  of 
the  constitutive  model  is  modified  to  include  tension  failure, 
shear  failure  and  an  effective  stress  law.  Tha  shape  of  the 
equivalent  source  spectrum  is  shown  to  be  strongly  dependent 
on  both  tension  failure  and  the  stress  relaxation  assumed  to 
result  from  the  effective  stress  formulation. 
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II.  DESCRIPTION  OF  MATERIAL  MODELS 

For  the  parameter  study  we  have  chosen  the  rock  con- 
stituent to  be  NTS  tuff.  The  following  sections,  2.1  and 
2.2,  present  the  pressure  and  material  strength  formulations 
used  in  the  parameter  study. 

2 . 1 THE  PRESSURE  COMPONENT 

Given  the  pressure  component  of  the  equation  of  state 
for  both  water  and  grain  density  tuff,  then  the  pressure  com- 
ponent for  a fully  saturated  mix  is  obtained  by  specifying 

M 

f * 5 — \--n-  * mass  fraction  of  water  (2.1) 

W M + M 

r w 

Figures  2.1  and  2.2  show  the  Hugoniot  and  release  isentropes 
for  these  two  constituents. 

At  a given  pressure  P a pressure  equilibrium  mix  is 
obtained  if 


v + f'v 
r w 

1 + f' 

e + f e 
r w 

1 + f' 


where 


(2.2) 


(2.3) 


(2.4) 


v , v , e and  e are  specific  volumes  and  energies  of  the 
rock  and  water  components  of  the  pressure  P.  It  is  now 
possible  to  generate  a table  that  gives  P as  a function  of 


I 


6 


Pressure  fiber) 


20 


0.2 


Figure  2. 


i iju : I— 

.22  0.24  0.3S 

Volume  (cn’/p) 

. Hugoniot  and  release  isentropes  for  tuff,  p ■ 

2.4  g/cc  [Riney,  et  aL  (1972)].  • 
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— f 21 

v and  e.  This  table  (Riney,  et  al . ) 1 J suitably  chooses 
values  of  v and  e so  that  no  time  consuming  search  is 
needed  when  performing  an  interpolation  at  a given  thermodyna- 
mic state.  Figure  2.3  shows  the  result  of  mixing  the  rock 
and  water  constituents  shown  in  Figs.  2.1  and  2.2  using  a 

0.17  value  for  f . 

w 

Air-filled  porosity  is  included  in  the  model  by  defining 
a parameter  a,  given  by 


V + v + V 
r w p m v 

V + V v 

r w v 


(2.5) 


where  v is  the  specific  volume  of  the  partially  saturated 
mix.  Then 


P - | P (v/a,  e) 


(2.6) 


where  P is  found  from  the  table  using  v,  e.  Equation 
(2.6)  is  equivalent  to  the  porosity  model  proposed  by 
Herrmann  (3)  except  P is  weighted  with  1/a.  Both  v and 
e are  zone  variables  and  are  obtained  from  SKIPPER  at  any 
given  time  cycle. 

The  parameter  a is  assumed  to  va.ry  in  the  following 
way  during  loading  (pore  collapse)  [Cherry,  ejb  al.]  ^ 


P > P 


a * 1 + (a  -1) 
e 


a ■ a + 6 
o 


‘P=0 


(2.7a) 

p-p  \ 

-Pc-U  P - Pc 

V 

(2.7b) 

n P/P  \ 

e ) 0 < P < P 
- - e 

(2.7c) 

< 0 (a'  = da/dp) 

(2.7d) 

9 


Pressure  (kba 
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where  P and  P are  the  elastic  and  crush  pressure.  B 
e c 

and  n are  determined  by  specifying  the  zero  pressure  bulk 

modulus  (k  ) of  the  porous  mix  and  requiring  that  a'  be 
o 

continuous  at  P - P • The  initial  air-filled  volume  fraction 

e 

(<{>  ) and  a are  related  by 
o o 


a - 1 


V + v + V 
w r p 


(2.8) 


while  the  volume  fractions  of  rock  and  water  (A  and  A) 

r w 

are  given  by 


1 + A 


- ( r-r^) 


(2.9) 


(2.10) 


where 


. r f 


Pr  ■ 2.4  g/cc 


Pw  * 1.0  g/cc 


We  have  found  (Cherry,  et  al^1  J ) that  a pore  recovery 
formulation  needed  to  be  included  in  the  model  in  order  to 
insure  a smooth  transition  between  loading  states  near  Pe» 
If  a*  is  the  minimum  a seen  by  a zone  then 


a = to*  + (1-t)  a (P) 


(2.11) 


where  (P)  is  given  by  Eqs.  (2.8)  and 


a"  - a 

t = 0,  a*  > a = , — — , 1 < a*  < 

— e l - a — — e 

ft 


(2.12) 


This  assures  full  recovery  when  a*  * ag  and  zero  recovery 


when  a*  ® 1. 
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Figure  2.4  shows  the  loading  and  release  P-v  curves 

that  result  from  adding  5 percent  air-filled  porosity  to 

the  mix  of  Fig.  2.3.  We  assume  P ■ 0.15  kbar  and  P = 

e c 

1.25  kbar.  Also  shown  for  comparison  in  the  figure  are  the 
hydrostatic  unloading  data  for  unit  3 tuff  at  the  Diamond 
Dust  site  (Stephens,  et  al . ^ ) . 

In  summary,  the  pressure  component  of  the  equation  of 
state  is  uniquely  determined  by  assuming  a pressure  equili- 
brium mix  and  specifying  f (Eq.  (2.1)),  <*>  (Eq.  (2.8)), 

o 

P , I , k (Eq.  (2.7b)  and  (2.7c))  and  the  grain  density  of 

6 C Q 

the  rock  component. 


2.2  MATERIAL  STRENGTH 

A strength  relation  of  the  form 


Sij  Sij  - 3 Y*  (PT'e) 


(2.13) 


where  S. . is  the  stress  deviator,  P = P + P , P equals 

xj  1 0 0 

the  overburden  pressure,  and 


Y(P,e) 


Y,  + Ym  F"  (2  ‘ r> 

m m 


U - ~)  P < Pm 
e m 

m 


(Y  + Y ) (1  - , P > P„ 

o n e — in 

m , 


(2.14) 


e > e 


m 


For  spherical  symmetry  Eq.  (2.13)  becomes 

|sj  £ f Y (PT,e)  (2.15) 

For  all  the  calculations  performed  for  the  parameter 
study,  a generalized  associated  flow  rule  was  used  for  the 
material  in  which  the  volumetric  plastic  strain  rate  was 
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assumed  to  be  zero.  This  assumption  implies  (Cameron  and 
Scorgie  ^ ) that  if  IsJ  exceeds  j Y,  then 

St  = (|  Y)  sign  (Sj)  (2.16) 

The  material  strength  dependence  on  P^ , , Y^  and 

P is  summarized  in  Fig.  2.5.  The  shock  loading  path  is 

in 

also  shown  in  the  figure  to  have  a slope  of  u/k  where  k is 
the  effective  bulk  modulus  of  the  shpck  wave.  This  loading 
path  is  only  valid  for  regions  where  the  radial  strain  is 
much  greater  than  the  tangential  strain,  i.e.,  at  the  shock 
front. 

In  Section  V the  material  strength  formulation  will  be 
modified  to  include  brittle  failure  and  an  effective  stress 
law. 


14 
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Yn  ■ Maximum  allowable  rtress  difference 

PT  = Hydrodynamic  component  oc  the  stress  tensor 
including  overburden  pressure 

« Overburden  pressure 

Figure  2.5.  Assumed  relationship  between  the  material  strength 
(Y)  and  the  hydrodynamic  component  of  stress  (PT) . 
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III.  DISCUSSION  OF  RESULTS 


3.1  THEORETICAL  CONSIDERATIONS 

A detailed  discussion  of  the  reduced  displacement 
potential/  RDP,  which  serves  as  the  equivalent  elastic  source 
representation  for  spherically  symmetric  explosions  is  pre- 
sented in  Appendix  A.  The  SKIPPER  code,  which  is  used  to 
propagate  the  shock  wave  into  the  elastic  region  is  discussed 
in  Appe-iix  B.  The  results  presented  in  Fig.  A.l  show  the 
calculated  reduced  velocity  potential  3pectra  using  the 
material  model  described  in  the  preceding  section.  As  noted 
in  Appendix  A,  the  spectra  are  essentially  flat  in  the  tele- 
seismic  frequency  band  for  device  yields  at  least  as  large 
as  200  kt.  This  result  along  with  Eq.  (A. 9)  imply  that  body 
wave  magnitude,  m^,  has  the  proportionality 

^ log  [cYM]  (3.1) 


where  («)  is  the  steady  state  value  of  the  RDP  and  c 
is  the  near  source  compressional  wave  speed.  Of  course, 
this  relation  does  not  include  the  effect  of  pP  or  any  of 
the  details  of  the  geophysical  properties  of  the  crustal 
layers.  It  has  been  derived  by  considering  the  impedenc* 
mixmatch  between  the  low  velocity  source  region  and  the 
basement  rock,  i.e.,  a t*.;o  layer  earth  model.  It  is  pre- 
fered  to  the  homogeneous  earth  assumption  which  gives 
'Moo)/o  for  proportionality. 

If  ^(w)  and  («)  are  two  values  of  ¥ («)  corres- 
ponding to  materials  i and  k then  the  change  in  tele- 
seismic  magnitude  is  g_von  by 


Am  ■ m 


i 


k 


m 


log 


(3.2) 
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I 


If  material  properties  remain  constant  and  the  only  change 
is  the  device  yield,  W,  then  cube  root  scaling  gives 

i k Wi 

Am  - in  - m*  - log  (3.3) 

Wk 


Equation  (3.2)  gives  the  "scaling"  law  for  teleseismic 
magnitudes  in  terms  of  cY(«).  This  equation  is  always  valid 
as  long  as  the  spectrum  of  the  equivalent  source  is  flat 
within  the  teleseismic  frequency  band  and  path  effects  asso- 
ciated with  earth  structure  are  invariant  between  events. 

It  should  be  noted  that  uncertainties  in  yield  esti- 
mates based  on  experimentally  determined  seismic  magnitudes 
may  be  related  to  uncertainties  in  near  source  material  pro- 
perties via  direct  application  of  Eq.  (3.2). 

3.2  PARAMETER  SENSITIVITY 

The  results  of  the  sensitivity  study  are  summarized  in 
Table  3.1.  Some  general  remarks  are  in  order  before  discus- 
sing the  effect  of  individual  parameter  variations.  Note 

that  the  ambient  density,  p , and  the  bulk  modulus,  k , are 

o o 

included  in  the  table  although  they  are  not  independent  para- 
meters in  the  context  of  this  study;  i.e.,  their  values  are 
determined  by  the  independent  parameters  appearing  in  the 
first  ten  columns.  Calculation  No.  18  is  selected  as  the 
reference  for  body  and  surface  wave  magnitude  variations  in 
order  to  avoid  the  appearance  of  minus  signs  in  the  Am^ 
column.  Calculations  marked  with  an  asterisk  in  the  Y (®) 
column  were  performed  with  a new  source  description  which 
gave  an  RDP  42  percent  higher  than  the  RDP  obtained  in 
Calculation  No.  8 using  the  old  source  description;  the  RDP's 
calculated  with  the  new  source  description  were  therefore 
scaled  down  by  a factor  of  0.704  and  the  scaled  values  are 
reported  in  Table  3.1.  All  calculations  appearing  in  Table 
3.1  were  performed  with  a device  yield  of  20  tons. 
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The  effects  of  water  content  on  seismic  coupling  are 
revealed  by  Calculations  1 through  5.  The  reduced  displace- 
ment potential  was  calculated  for  five  fully  saturated  tuffs 
having  different  mass  fractions  of  water,  f . Fully  saturated 
tuffs  were  chosen  in  order  to  eliminate  the  effects  of  air- 
filled  porosity.  Note  that  adding  water  decreases  both  the 

density,  p , and  the  bulk  modulus  of  the  mixture,  k . Since 
o o 

the  shear  modulus  is  held  constant  at  40  kbars,  variations 

in  compressional  wave  speed,  c,  reflect  only  the  changes  in 
density  and  bulk  modulus.  Variations  in  body  wave  magnitude 
as  functions  of  water  mass  fraction  are  described  in  Fig. 

3.1.  The  decrease  in  seismic  coupling  as  the  water  content 
is  increased  is  to  be  expected.  Since  water  is  more  compres- 
sible than  rock  at  all  pressures  of  interest,  the  amount  of 
energy  dissipated  in  the  form  of  heat  via  shock  loading 
increases  with  increasing  water  content,  leaving  less  energy 
available  for  conversion  into  seismic  radiation,  m^  de- 
creases by  about  0.5  magnitude  units  as  the  mass  fraction 
of  water  increases  from  0 to  0.19. 

The  effect  of  air-filled  voids  on  seismic  magnitude 
is  obtained  from  Calculations  4,  6 and  7.  The  water  content 
is  held  constant  at  fw  * 0.17  for  these  and  the  remaining 
calculations  in  Table  3.1.  The  variation  of  seismic  magni- 
tude as  a function  of  the  void  fraction,  4)  , is  presented 

0 

in  Fig.  3.2.  Note  that  the  intrpduction  of  air-filled 
voids  causes  a drastic  reduction  in  seismic  coupling  and  the 
importance  of  establishing  this  parameter  preshot  is  obvious. 

Calculations  8 and  9 show  that  seismic  magnitude  is  not 

sensitive  to  P , the  pre^ jure  below  which  pore  closure  is 
© 

elastic  and  fully  recoverable.  A factor  of  two  increase  in 

P increases  the  seismic  magnitude  by  0.04  units.  Similarly, 
e 

Calculations  10  and  11  show  that  magnitude  is  insensitive  to 

variations  in  P , the  pressure  at  which  all  ambi'T.c  porosity 

c 

is  irreversible  removed.  An  increase  in  P„  by  a factor  of 

c 


1 
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2.5  produces  an  increase  in  magn.tude  of  0.03  units.  It 

should  be  noted  that  and  are  associated  with  the  porosity 

model  and  might  increase  ia  importance  as  <t>  , the  air-filled 

o 

void  fraction,  is  increased. 

The  effect  of  variations  in  near  source  seismic  veloc- 
ity, c,  on  magnitude  is  obtained  from  Calculations  8,  12  and 
13.  The  permissible  range  of  c is  limited  at  this  point 
because  the  properties  of  tuff,  the  properties  of  water,  the 
mass  fraction  of  water  and  the  shear  modulus  are  held  con- 
stant. Variations  in  c therefore  reflect  only  a change  in 
the  initial  slope  of  the  crush  curve  associated  with  the 
porosity  model.  As  shown  in  3,  a 6 percent  error  in 

P wavr  velocity  produces  a cha  m^  of  about  0.05  magni- 

tude units.  Increasing  the  compressional  wave  speed  increases 
seismic  coupling  for  body  waves. 

Calculations  8,  14  and  15  illustrate  the  effect  of 
variations  in  the  shear  modulus  on  seismic  magnitude.  Body 
wave  magnitude  is  reduced  as  the  shear  modulus  increases  as 
shown  in  Fig.  3.4. 

The  yield  surface  is  described  by  the  parameters  Y , 

Ym,  and  Pm.  Calculations  8,  16,  17  and  18  show  the  effect 
°f  Y , the  zero  pressure  yield  strength,  on  seismic  coupling. 
The  variation  of  Cm  with  Y^  is  plotted  in  Fig.  3.5.  In- 
creases in  Y^  tend  to  reduce  seismic  coupling.  The  results 
of  Calculations  7,  8 and  22  are  plotted  ir.  Fig.  3.6  which 
shows  the  effect  of  Yj^^,  the  maximum  allowable  material 
strength,  on  seismic  magnitude.  Increasing  Y by  a fac- 
tor  of  two  causes  magnitude  to  decrease  by  0.192  magnitude 
units,  indicating  that  Y^^  is  a fairly  sensitive  parameter. 

For  a given  value  of  Y^^,  the  shape  of  the  material 

strength  surface  is  controlled  by  P , the  pressure  at  which 

m 

YMAX  atta^ned-  From  Calculations  10  and  19  we  see  that 
increasing  P^  by  a factor  of  2.08  causes  magnitude  to  in- 
crease by  0.26  units.  This  result  indicates  that  the  shape 
of  the  failure  surface  plays  an  important  role  in  determining 


. 
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coupling  efficiency.  The  material  has  less  strength  at  low 
pressures  when  Pm  is  increased  and  we  should  expect  a cor- 
responding increase  in  magnitude.  In  fact,  if  any  of  the 

parameters  Y , Y , or  are  varied  such  that  the  material 
o m m 

strength  is  enhanced,  a decrease  in  coupling  efficiency  is 
observed . 

The  effect  of  overburden  pressure,  P , on  magnitude, 

o 

as  determined  by  Calculations  6,  10,  20  and  21,  is  shown  in 
Fig.  3.7.  If  all  material  vroperties  are  held  constant, 
coupling  efficiency  decreases  with  increasing  overburden 
pressure.  An  increase  in  the  depth  of  burial,  DOB,  by  a 
factor  of  six  causes  magnitude  to  decrease  by  a factor  of 
two.  The  decrease  in  magnitude  is  not  uniform  with  depth. 

An  increase  in  DOB  at  shallow  depths  causes  magnitude  to  de- 
crease more  than  that  caused  by  an  equivalent  increase  in 
DOB  at  greater  depths.  This  effect  results  from  the  rapid 
increase  in  material  strength  at  low  pressures.  We  should 
note  that  it  is  unlikely  that  material  properties  would  re- 
main invariant  as  overburden  pressure  is  increased  from  100 
bars  to  600  bars.  An  increased  DOB  usually  corresponds  to 
an  increase  in  water  content,  resulting  in  decreased  air- 
filled  porosity  and  material  strength. 
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Figure  3.7.  Effect  of  overburden  pressure  on  seismic  magnitude 
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IV.  SUMMARY  OF  RESULTS 


The  results  of  the  parameter  stucy  provide  valuable  in- 
sight concerning  the  importance  of  various  material  parameters 
with  respect  to  accurate  prediction  of  seismic  coupling.  It 
has  been  shown  that  increasing  the  overburden  pressure,  or, 
equivalently,  the  depth  of  burial,  decreases  coupling  effi- 
ciency substantially.  Therefore,  the  DOB  muse  be  accurately 
known  in  order  to  make  a reliable  prediction  of  seismic  magni- 
tude. It  is  also  desirable  to  know  the  properties  of  the 
material  in  the  immediate  vicinity  of  the  device;  if  sub- 
stantial spatial  variation  of  the  material  properties  occurs 
near  the  device  or  if  the  DOB  is  quite  shallow  it  may  be 
necessary  to  consider  two-dimensional  effects  which  are  not 
treated  here.  Accurate  prediction  of  seismic  magnitudes 
via  one-dimensional  spherically  symmetric  calcul«itions  is 
therefore  understood  to  be  dependent  on  the  material  uni- 
formity about  the  device. 

Since  the  elastic  properties  of  the  near  source  mate- 
rial appear  in  the  magnitude  relation,  Eq.  (3.2),  via  the  P 
wave  velocity,  c,  it  is  obviously  important  to  determine 
their  values  accurately.  A positive  error  in  determination 
of  c would  cause  body  wave  magnitude  to  be  over-predicted. 

A positive  error  in  y would  cause  magnitude  to  be  under- 
predicted. 

The  steady  state  value  of  the  RDP,  ¥(»),  which  appears 

in  the  magnitude  relation,  Eq.  (3.2),  is  dependent  on  the 

shock  response  of  the  near  source  material.  The  near  source 

material  properties  which  have  the  most  pronounced  effect 

on  the  shock  response  are  the  water  macs  fraction,  fw,  and 

the  air-filled  void  fraction,  <p  . lositive  errors  in  either 

s 

f or  # could  lead  to  substantial  under-prediction  of 

w g 

seismic  magnitudes.  Seismic  magnitude  is  not  very  sensitive 
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to  P , the  elastic  pressure#  and  P , the  crush  pressure, 

6 C 

indicating  that  details  in  the  porosity  model  are  relatively 
unimportant. 

Seismic  magnitudes  are,  however,  very  sensitive  to  the 

parameters  which  describe  the  failure  surface.  If  Y#,  Ym, 

or  P are  varied  such  that  the  material  strength  is  enhanced 
m 

the  coupling  efficiency  of  an  explosive  device  is  impaired. 
Thus,  a positive  error  in  the  material  strength  would  lead 
one  to  an  under-prediccion  seismic  magnitude. 
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V.  THE  EFFECT  OF  MATERIAL  STRENGTH  ON  THE  PRESSURE  AND 
DEVIATORIC  COMPONENTS  OF  THE  STRESS  TENSOR 

5.1  INTRODUCTION 

Most  code  calculations  of  the  type  presented  in  Sec- 
tion 3,  either  for  ground  shock  predictions  or  equivalent  source 
calculations,  have  assumed  that  material  strength  depends  only 
on  the  pressure  component  of  the  stress  tensor  (i.e.,  Eq.  (2.14)). 
Also,  very  little  attention  has  been  paid  to  the  consequences 
of  exceeding  the  allowable  material  strength  other  than  envok- 
ing  either  an  associated  or  non-associated  flow  rule,  i.e., 

Eq.  (2.16).  Further,  the  influence  of  water  content  on  mate- 
rial strength  via  an  effective  stress  law  has  not  been  considered 
in  any  ground  shock  calculation  performed  to  date. 

A number  of  investigators  have  considered  the  above 
topics  separately.  In  particular.  Cherry  and  Petersen^ 
have  identified  the  dependence  of  material  strength  on  the 
third  deviator ic  invariant.  Maenchen  and  Sack^10^  presented 
a formulation  for  tension  failure  which  is  ideally  suited 
for  a stress  wave  code.  Garg  and  Nur^^  investigated  the 
consequences  of  various  effective  stress  laws  on  the  inter- 
pretation of  experimental  data.  We  have  relied  heavily  on 
the  results  of  these  investigators  in  order  to  formulate  the 
strength  portion  of  the  constitutive  model. 

5.2  THE  DEPENDENCE  OF  MATERIAL  STRENGTH  ON  THE  STRESS  IN- 
VARIANTS  

Figure  5.1  is  taken  from  Cherry  and  Petersen^  and 

shows  the  improvement  in  the  definition  of  material  strength, 

Y,  when  P is  replaced  by  P.  Similar  improvement  has  also 

been  found  for  granite,  limestone  and  glass.  Y,  P and  P are 

obtained  from  the  stress  invariants  J , J' , J",  as  follows: 

1 2 ) 
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Figure  5.1.  These  figures  are  from  Cherry  and  Petersen  [1970] 
and  show  the  improvement  in  failure  surface  de- 
finition when  P is  replaced  by  P.  The  strength 
data  are  from  Handin,  et  al.  [1967]. 
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Y = ( 3J  ) 
2 


p ■ p 


- S fr) 


J'\  1/2 


(5.2) 


(5.3) 


P is  now  the  pressure  component  of  the  stress  tensor  discussed 
in  Section  II  with  the  overburden  pressure,  P^,  added. 

J'  and  J'  are  the  first,  second  deviatoric  and  third  devia- 

2 i 

toric  stress  invariants.  If  o , o , o are  principal 

11  22  11 

stresses,  then 


J * o +o  +o 
1 11  22  11 


(5.4) 


(a  +P) 2 + (c  +P)  2 + (o  +P) 2 


(5.5) 


J'  ■ (a  +P)  (o  +P)  (o  +P) 
1 11  22  11 


(5.6) 


Note  that  when  the  intermediate  principal  stress,  o^,  is 
equal  to  either  the  maximum,  a , or  minimum,  c , principal 

ii  ii 

stress, then 


Y - la  - a 

11  2 2 


p . . (°»»  i g») 


(5.7) 


(5.8) 


The  functional  form  we  have  chosen  to  represent  the 
variation  of  material  strength,  Y,  with  the  stress  state,  P, 
is  identical  to  that  used  in  Section  2,  i.e., 
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(Y  + Y ) 
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P > P 


m 


(5.9) 


e < e 


m 


Y - 0 


e > e 


rn 


where  Y , Y , P and  e are  input  constants.  These  con- 
o m m m 

stants  along  with  P , the  overburden  pressure,  determine 

0 

the  variation  of  material  strength  with  stress  state  along 
with  the  initial  stress  state,  P , corresponding  to  the 
overburden  pressure  at  the  depth  of  burial  of  the  device. 


The  material  strength  dependence  on  P , Y , Y and 
_ e o m 

Pm  is  summarized  in  Fig.  5.2.  The  shock  loading  path  is 

also  shown  in  the  figure  to  have  a slope  of  l-2c  , where  a 

is  Poisson's  ratio.  This  loading  path  is  only  valid  for 

regions  where  the  radial  strain  is  much  greater  than  the 

tangential  strain,  i.e.,  at  the  shock  front. 


5.3  SHEAR  FAILURE 

Hooke's  law  is  used  to  obtain  an  initial  estimate  of 
the  stress  deviators,  i.e., 

sjj1  - + 2u  eij  At  (5.10) 

where  and  are  the  values  of  the  stress  deviator 

at  time  t + At  and  t respectively. 

u is  the  shear  modulus 

e^j  is  the  strain  rate  devaitor 

At  is  the  time  increment 


| i 

i 


i 


i 
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Shear  failure  occurs  if  the  material  strength,  given 
by  Eq.  (5.9)  and  evaluated  at  P,  is  exceeded,  i.e.,  if 


J > Y (P) 


(5.11) 


where 


P - P 
,n+l 


| (s2  + S2  + Ss  )] 1/2 

2 \ 11  2 2 S S 'J 

(A  A A \ 

8,1  "»  S 


(5.12) 


n+1 


1/3 


(5.13) 


:n+i 


P includes  the  overburden  pressure  and  we  have  omit- 

A 

n+1  superscript  on  S^. 

If  J > Y (P)  then  adjustment  of  the  stress  deviators. 


, is  required.  We  assume  that 


Sn+^"  * a Sn+^ 
sij  a sij 


(5.14) 


,n+l 


where  is  the  adjusted  value  of  the  stress  deviator 

and 


a 


Y(P)  + | 

(:,  *,  >,r 

\ 2 / 

j + 1 ( 

s s s \ 1/3 

1 2 I 

2 

(5.15) 


dY 


b » — (evaluated  at  P) 
dP 


(5.16) 


Equation  (5.15)  is  obtained  by  approximating  the  strength 


function  at  P with  a first  order  Taylor  series  and  assum- 


ing that  no  adjustment  in  pn+^  occurs  during  shear  failure. 


I 
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5.4 


TENSION  FAILURE 


Tension  failure  is  assumed  to  occur  in  the  element  if 
a principal  stress  is  greater  than  zero  and  if  the  inequality 
of  (5.11)  has  evei  been  satisfied.  We  then  apply  the  tension 
failure  model  proposed  by  Maenchen  and  Sack*’10*  and  introduce 
an  inelastic  strain  normal  to  the  crack.  This  inelastic  strain 
is  just  sufficient  to  zero  the  tensile  stress. 


For  example,  if  c , o and  a are  the  three  princi- 

/v  11  2 2 ) S 

pal  stresses  and  if  o is  greater  than  zero,  then  the  ad- 


1 1 


justed  stress  (c  , a , a ) are  given  by 

11  22  ss  1 


o ■ o - (k  + i y)  ie 
ii  li  3 ii 


22 


o - (k  - -r  y)  Ae 
2 2 3 1 1 


(5.17) 


3 3 


o - (k  - T u)  Ae 

3 3 3 n 


where 


Ae 


1 1 


k + j u 


(5.18) 


k is  the  bulk  modulus,  u is  the  shear  modulus  and  all 
stresses  include  the  overburden  pressure. 

If  two  of  the  principal  stresses,  say  a and  a 

112  2 

are  greater  than  zero,  then  the  stress  adjustment  becomes 


o =o  - (k  - T u)  (Ae  + Ae  ) - 2|iAe 
11  11  3 11  22  11 


o = o - (k  - T u)(Ae  + Ae  ) - 2uAe  (5.19) 

22  22  3 11  22  22 


o =0  - (k  - u)  (Ae  + Ae  ) 

33  S3  3 11  22 
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where 


Ae 


1 1 


Ae 


2 2 


(k  + ^ ) a — (k  - 3-  u)  a 

£ U £ LL 

4y  (k  + J) 

(k  + ^ u)  o - (k  - - y)  o 

£ LI £ LL 

4y  (k  + 


(5.20) 


All  the  inelastic  strain  increments  (Ae  , Ae  , Ae  ) are 

11  22  33 

accumulated  on  each  cycle  giving 


,n+I  _r,  . * 

: = E + Ae 

li  ii  li 


„n+l  _n  , A 
l *=  E + Ae 
22  22  22 


(5.21) 


_n+l  _n  ^ . 

E = E + Ae 

33  33  33 

These  equations  give  the  basic  stress-strain  adjustment 
during  tension  failure.  However,  they  apply  equally  well  for 
crack  closure.  If  at  least  one  of  the  total  strains  in  Eq. 
(5.21)  (say  En  ) is  greater  than  zero,  then  the  crack  will 

ii  * * 

open  or  close  depending  on  the  sign  of  a . If  a >0, 
then 


Ae  >0 
1 1 


En+1  > e" 
11  11 


and  the  crack  width  increases.  The  inequalities  are  reserved 

A 

if  T <0  and  the  crack  width  decreases.  Closure  will  con- 
1 1 

tinue  until 


En  + Ae  < 0 . 

ii  ii 
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Then 


and  the  crack  is  completely  closed.  When  this  state  is 
achieved  the  element  is  able  to  support  a compressive  stress 
in  the  (1,0,0)  direction. 
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5.5 


EFFECTIVE  STRESS 


All  effective  stress  laws,  formulated  for  fluid- 
saturated  porous  rocks,  assume  that  the  fluid  pressure  should 
be  included  in  the  calculation  of  the  stress  state  which  de- 
termines material  strength.  An  excellent  review  of  various 
definitions  cf  effective  stress  and  their  consequences  in 
terms  of  laboratory  strength  data  has  been  given  by  Garg  and 
Nur. [11! 

The  definition  of  stress  state  we  have  adopted  is  P, 
where  P is  given  by  Eq.  (5.2).  If  Pp  is  the  pore  fluid 
pressure,  then  the  effective  stress  state  becomes 

J'  ■ 1/3 


- p(pp) 


where  F(Pp)  is  a positive  function  of  the  pore  fluid  pres- 
sure. We  have  assumed  that  F(Pp)  behaves  as  follows: 


F(P  ) = 0 if  tension  cracks  are  open. 
P 


(5.22) 


F(p  ) « 0 if  P has  never  exceeded  the  crush 
P pressure,  Pc» 


(5.23) 


r * v 1/3 


F (P  ) = 0 if  P 
P 


-i(?)  <o- 


(5.24) 


If  none  of  the  above  apply,  then 


F(Pp) 


-*(?) 


1/3 


(5.25) 


The  conditions  given  in  Eqs.  (5.22)  and  (5.23)  assume 
that  the  pore  fluid  pressure  has  no  effect  on  the  stress  state 
when  air  filled  porosity  exists  in  the  element.  If  no  air  filled 
porosity  exists,  then  Eqs.  (5.24)  and  (5.25)  apply. 
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Figure  5.3  shows  the  effect  on  material  strength  when 
P > Pc.  The  relaxation  to  Y#  is  achieved  by  assuming  the 
material  behaves  like  a Maxwell  solid.  The  stress  deviator 
S,  . is  given  by 


an+1 

Sij 


H1- 


J > Y 


J < Y 


where  S?+1  and  J are  determined  using  Eqs.  (5.10)  and 

(5.12).  The  relaxation  time,  6/(l-Y  /J) , increases  as  J 

e 

approaches  Y . This  gives  a smooth  transition  between  the 
Maxwell  solid  and  elastic  behavior. 

The  behavior  of  F(P  ) and  the  subsequent  relaxation 
process  given  by  Eqs.  (5.22)  through  (5.26)  is  certainly  an 
oversimplification  of  the  role  of  pore  fluid  pressure  in 
determining  effective  stress.  Future  work  o i this  topic 
should  merge  experimental  data  of  the  type  obtained  by  Duba, 
et  al. with  the  theoretical  work  of  Garg  and  Nur. 

We  do  feel,  hv-vever,  that  the  main  features  of  effective  stress 
are  presently  contained  in  the  model.  The  primary  input  para- 
meters t<hich  control  the  process  are  Y , P„  and  6.  At 

6 CT 

least  some  ipsight  concerning  the  appropriate  values  for  Y# 
and  P ' should  be  able  to  be  obtained  from  laboratory  experi- 


ments. The  value  assumed  for 


controls  the  relaxation 


time  and  we  have  used  close-in  tine  history  data  to  determine 
this  parameter. 


f 


Figure  5.3.  The  relaxation  of  material  strength  when  all 
air  filled  porosity  has  been  removed. 
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5.6  SUMMARY  OF  INPUT  PARAMETERS  FOR  MATERIAL  STRENGTH 

The  limiting  value  of  material  strength,  Y (f>) , is 
obtained  by  specifying  Y , Y . F and  e^.  The  function  Y (F) 
is  given  in  Eq.  (5.11).  Except  for  the  definition  of  P, 
these  parameters  and  Y (F)  are  identical  to  that  of  Section 
2.2.  The  overburden  pressure  P#  is  required  to  determine 
the  ambient  stress  state  of  the  material  which  is  also  re- 
quired in  Section  2.2.  The  crush  pressure,  P , plays  a dual 
role.  First  it  determines  the  limiting  pressure  at  which 
the  steepest  release  path  occurs  in  the  P-V  plane,  as  shown 
in  Fig.  2.4.  Second,  it  establishes  the  pressure  at  which 
an  effective  stress  law  is  used  to  determine  material  strength 
as  shown  in  Fig.  5.3.  The  parameter  6 determines  the  re- 
laxation time  of  the  stress  deviators  during  the  application 
of  the  effective  stress  law. 

It  is  worth  emphasizing  that  the  only  new  parameter 
introduced  in  this  section  is  6,  the  relaxation  time.  How- 
ever, the  interpretation  of  material  failure  has  changed 
considerably  over  the  simple  plastic  flow  model  used  in 
Section  2,  Eq.  (2.16).  We  now  present  rome  of  the  conse- 
quences of  this  revision. 
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5.7  THE  CONSEQUENCES  OF  MATERIAL  FAILURE  ON  THE  EQUIVALENT 

ELASTIC  SOURCE 

A number  of  equivalent  source  calculations  were  per- 
formed using  the  material  strength  constitutive  model  pre- 
sented in  Sections  5.2  through  5.5.  Table  5.1  lists  the  key 
material  parameters  used  in  the  calculations  along  with  the 

calculated  values  of  ¥ (®)  and  the  final  cavity  radius,  R , 

1 cav 

for  a device  yield  of  0.02  kt. 

The  input  data  listed  in  Table  5.1  has  been  taken  from 
a number  of  sources.  In  particular  the  rock  mechanics  labora- 
tory data  of  Stephens,  Schock  and  Heard  at  the  Lawrence  Liver- 
more Laboratory  were  used  for  granite  and  rhyolite . The  data 
taken  by  Terra  Tek  for  many  events  at  NTS  Area  12  were  used 
for  tuff.  We  should  note,  however,  that  no  material  strength 
data  could  be  found  for  saturated  rhyolite.  We  assumed  that 
the  material  strength  of  this  rock  type  would  be  the  same 
as  granite. 

Figures  5.4  and  5.5  show  velocity  time  histories  re- 
corded at  two  distances  from  Piledriver,  a 60  kt  event  in 
granite  at  NTS  (Parret^^).  In  the  past,  most  calculations 
using  measured  granite  strength  data  have  been  unsuccesful 
in  duplicating  these  velocity  histories. 

Figures  5.6  and  5.7  show  the  velocity  histories  ob- 
tained from  our  calculation  of  the  Piledriver  event  (Calcula- 
tion 130)  and  Table  5.2  summarizes  the  comparison  between  the 
calculation  and  the  data.  Cavity  radius  and  peak  velocities 
at  the  two  stations  agree  with  the  data  to  within  5 percent. 
However,  residual  displacement  at  our  first  station  is  too 
high.  This  is  due  largely  to  the  difficulty  in  reproducing 
the  large  negative  phase  of  the  velocity  pulse  of  Fig.  5.4. 
Since  this  difficulty  has  been  encountered  ii>  all  previous 
calculations  of  Piledriver,  it  was  decided  to  concentrate 
on  the  second  station.  Here  we  attained  excellent  agreement 
with  the  velocity  data,  including  the  long  duration  negative 
phase  of  the  pulse.  This  gave  U3  good  agreement  with  the 
residual  displacement  measurement  as  well. 
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TABLE  5.2 

COMPARISON  BETWEEN  DATA  AND  CALCULATION 
FOR  PILEDRIVER 


Cavity  radius  (meters) 

Peak  Velocity  ft/sec 
Station  1 (660'  range) 
Station  2 (1543'  range) 

Residual  displacement  (inches) 
Station  1 
Station  2 


Data  Calculation 
44.5  46.5 

110  103 

18.8  19.4 

22  38 

7.2  6.1 


Velocity  history  at  668  feet  horizontal  range  from  Piledriver 


Figure  5.5.  Velocity  history  at  1543  feet  horizontal  range  from  Piledirve 
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Figure  5.6.  Calculated  particle  velocity  at  668  feet  from 
Piledriver . 
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Figure  5.7.  Calculated  particle  velocity  at  1548  feet 
from  Piledriver. 
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The  Hoggar  Granite  calculation  was  intended  to  be  re- 
presentative of  the  source  coupling  for  the  French  nuclear 

tests  conducted  in  the  Sahara.  For  this  environment  P was 

c 

assumed  infinite,  a value  which  corresponds  to  a dry  material. 
From  Table  5.1  we  find  that  the  calculated  cavity  radius  for 
Piledriver  was  3.2  meters  while  for  the  Hoggar  granite  the 
cavity  radius  was  1.9  meters.  This  gives  a factor  of  five 
difference  in  cavity  volumes  between  the  two  test  areas  which 
is  consistent  with  the  published  French  test  results. 

Figure  5.8  shows  the  source  functions  corresponding 

to  Calculations  127,  130  and  125.  Both  the  Piledriver  and 

Pahute  Mesa  spectra  are  noticeably  peaked.  The  amplitude 

of  the  peak  seems  to  correlate  directly  with  Y . 

o 

In  order  to  identify  the  physical  phenomenon  responsible 
for  the  peak  we  show  the  two  Hoggar  source  spectra  correspond- 
ing to  Calculations  128  and  129  in  Figs.  5.9  and  5.10.  For 
Calculation  128  neither  stress  relaxation  or  tension  failure 
was  allowed  and  the  source  spectrum  is  flat;  a result  which 
is  identical  to  the  spectra  shown  in  Fig.  A. 2 in  Appendix  A. 
Tension  failure  was  permitted  in  Calculation  129  and  the 
spectral  peak  for  this  calculation  is  approximately  a factor 
of  three  greater  than  Y {») . 

The  effect  of  stress  relaxation  on  spectral  shape  is 
obtained  by  comparing  the  source  spectra  from  Calculations 

124  and  125.  The  total  relaxation  is  less  for  124  than  for 

125  since  the  corresponding  values  of  P are  0.5  kbar  and  1.0 

c 

kbar.  Figure  5.11  compares  the  source  spectra  from  these  two 
calculations  and  indicates  that  the  greater  the  stress  relaxa- 
tion the  larger  the  spectral  peak  when  compared  to  V (®) . 

Both  tension  failure  and  stress  relaxation  within  the 
inelastic  region  seem  to  be  the  primary  physical  mechanisms 
responsible  for  producing  a peaked  spectrum.  There  is  a 
good  explanation  for  this  since  both  mechanisms  introduce  a 
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Figure  5.9.  Equivalent  source  function  for  Calculation  128. 


R-2742 


discontinuity  in  the  magnitude  of  allowable  stress  between  the 
nonlinear  and  linear  elastic  regions.  This  discontinuity 
causes  the  velocity  field  in  the  elastic  region  to  rarefy, 
i.e.,  to  become  negative.  The  negative  portion  of  the 
velocity  time  history,  shown  in  Fig.  5.7,  causes  the  peak  in 
the  source  spectrum. 
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Figure  5.11.  Equivalent  source  spectra  for  Calculations 
124  and  125. 
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APPENDIX  A 


THEORETICAL  SCALING  OF  m.  WITH  REDUCED 

D 

DISPLACEMENT  POTENTIAL 


The  equivalent  elastic  source  for  a spherically  sym- 
metric explosion  is  often  represented  by  a reduced  displace- 
ment potential,  4'(t).  This  quantity  completely  characterises 
the  coupling  of  the  explosion  into  elastic  waves  in  the  im- 
mediate source  region.  Therefore,  we  expect  proportionality 
between  ¥ and  teleseismic  waves.  While  this  proportionality 
can  be  quite  complicated,  it  is  useful  to  have  a simple 
scaling  law  which  indicates  the  variation  of  teleseismic  body 
wave  amplitude  with  4*  and  the  properties  of  the  emplace- 
ment material.  Such  a scaling  relationship  is  derived  in 
this  appendix. 

We  begin  by  reviewing  the  most  elementary  model,  that 
of  a delta  function  source  in  a one  material  or  homogeneous 
earth  model.  We  then  add  material  inhomogeneity  in  the  form 
of  a second  material  and  allow  a more  complex  source  repre- 
sentation, arriving  at  the  scaling 

|utl  * V-  (A‘] 

where  |u  | is  the  amplitude  of  some  consistently  measured 

# 0 

phase  on  the  teleseismic  record  and  is  a constant  which 

characterizes  the  amplitude  of  the  explosion  generated  pulse. 


A Spherically 


letric  Source  in  a Homogeneous  Earth 


A spherically  symmetric  point  source  of  dilatational 
waves  can  be  characterized  by  a reduced  displacement  poten- 
tial. That  is,  the  radial  displacement  u (R,t)  can  be 

s 

written 


ug (R, t ) 


41  (t  ) 
Ra 
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whire  t * t - R/ag,  the  time  retarded  by  the  P wave  travel 
time  to  the  position  R.  Here  ag  is  the  P wave  velocity  in 
the  source  region. 

Taking  Fourier  transforms,  (A. 2)  may  be  written 


ug 


fc  * *) 


¥(u)  e 


-ik  R 
a 


(A. 3) 


where  ka  = w/ag.  If  kQR  >>  1,  that  is,  if  R is  many  wave- 
lengths, only  the  "far  field"  terms  need  to  be  retained  and 


ug  (R, u>) 


V(m)  "ikaR 
Ras  B 


(A. 4) 


For  teleseismic  body  wave  amplitudes  only  frequencies 

less  than  about  2.5  Hz  are  of  importance.  Then  for  low  yield 

explosions  ¥ (w)  is  essentially  constant  and  equal  to  V , 

* " 

the  static  value  of  the  reduced  displacement  potential. 
Therefore,  for  a low  yield  source  of  the  type  discussed  in 
Section  IV,  and  a homogeneous  earth  model,  we  have 


- la 

% ,T 


(A. 5) 


For  higher  yields  the  teleseismic  frequency  band  in- 
cludes portions  of  the  peak  on  the  4* (u>)  curves.  Therefore 
Eq.  (A. 5)  is  improved  by  replacing  4^  by  an  "equivalent" 
value,  4^,  which  is  yield  dependent.  Once  again,  this  is 


The  assumption  that  4'<w)  = 4,oo  for  the  teleseismic  frequen- 
cies is  equivalent  to  assuming  that  the  source  appears  to  be 


a delta  function  of  amplitude 


to  a teleseismic  observer. 
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equivalent  to  assuming  that  the  source  behaves,  teleseis- 
mically,  as  if  it  were  a delta  function  of  amplitude  4'*. 
Computation  of  theoretical  seismqgrams  show  that  this  assump- 
tion's not  bad.  The  value  of  4e  is  taken  to  be  the  val-ie 
of  ¥ (w)  at  the  frequency  which  cube  root  scales  to  the 
yield  of  interest.  With  this  modification,  (A. 5)  is  the 
appropriate  scaling  law  for  the  sources  of  Section  IV  in  a 
homogeneous  earth  model. 

A Spherically  Symmetric  Source  in  a Two  Layered  Whole  Space 


Let  us  now  assume  that  the  spherically  symmetric  source 
is  located  in  a space  composed  of  two  layers.  The  geometry 
is  shown  in  Fig.  A.l. 


^Source 


s (a#,  p^) 


‘V  V °b) 


Figure  A.l.  A source  in  a two  layered  medium. 


If  attention  is  restricted  to  the  far  field,  the 
displacement  in  potential  form  is  given  by  (A. 4).  Then  the 
P wave  displacement  at  the  point  B below  the  interface  is 
given  by 

- A ¥(w)  ”iwTb 
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where 


(v * is)  (V * £■») 


2(5  a / \ 

a - -F p,  (8bf,x  + Bsp.Y)  - 

o ■ °«ab6s6b  * W.V’  + ‘,s6sp,^Y, 
c 

♦ >Vb  (VbV.  + as6bp,p.)  + ^7  ViV- 


X » P, 


- Y ’ #. + ?r  1 ” °b " ‘ fr 


(A. 7) 


p ■ COS0 
1 


/ 8*\l/2  / f>v  \1/2 

,'p,-\1-^)  ' p,  - C0SV  p.  * i1  - ct)  •' 


8 „ b 

c “ sin0s  sin0fa 


. q - 2 (obeJ  - 0,6*) 


In  (A. 6)  the  exponential  tern  is  merely  the  phase  shi 

due  to  the  travel  tine  iron  the  source  to  B.  The  L is  t e 

, wave.  For  a homogeneous 

spreading  tern  for  the  spherical  wave,  to 

material,  L - A,  the  distance  between  the  source  and  B.  The 
fficient  A is  the  refraction  coefficient  which  depen  s 
Sertiel  of  the  two  media  and  the  angle  of  incidence 

which  we  denote  0g* 

Application  to  Teleseismic  Body  Wave_s 

The  two  layered  model  of  Fig.  A.l  is  more  realistic 
than  a homogeneous  earth  model  for  illuminating  the  gross 
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properties  of  the  teleseismic  body  waves  excited  by  explosions 
in  different  source  materials.  For  example,  consider  the  NTS 
explosions.  The  resulting  teleseismic  body  waves  observed  at 
a single  seismograph  station  are  expected  to  vary  due  to  dif- 
ferences only  in  the  explosion  source  and  near  source  crustal 
structure.  All  other  features  of  the  travel  path  are  essen- 
tially held  constant. 


Below  some  depth  one  expects  all  of  the  NTS  to  be 
underlain  by  a common  material.  The  material  denoted  by  a 
subscript  b can  be  thought  of  as  this  basement  rock.  Ex- 
plosions at  NTS  have  been  detonated  in  such  widely  varying 
materials  as  alluvium,  tuff,  rhyolite  and  granite.  These  are 
characterized  by  P wave  velocities  ranging  from  2.0  - 5.5 
km/sec.  The  subscript  s can  be  taken  to  represent  this  source 
material.  For  the  long  wavelengths  (>  2 km)  of  interest  for 
teleseismic  waves,  the  single  interface  transition  between 
the  source  and  basement  materials  is  a useful  approximation. 

The  body  waves  reaching  the  teleseismic  field  exit  the 
source  region  at  near  vertical  takeoff  angles.  That  is, 
cos6g  z 1 and  cosO^  » 1.  Then  it  is  clear  from  (A. 7)  that 


L 


a. 


(A. 8) 


However,  since  th~  basement  rock  is  held  in  common  from  event 
to  event,  one  can  choose  arbitrarily  large  while  Rg  re- 
mains rather  small.  Thus,  comparing  different  events  at  NTS, 

L s a, /a  . 
o s 

The  complexity  of  its  form  makes  it  difficult  to  de- 
duce a similarly  simple  relationship  for  A.  However,  numeri- 
cal experiments  using  values  for  typical  geologic  materials 
and  horizontal  wave  speed  c (or  takeoff  angle  0^)  show  that  A 
is  ”sry  nearly  proportional  to  o /a.  . Taking  c = 12  km/sec 

9 O 

and  ■ 6.0  km/sec,  *=  3.5  kn/sec,  = 2.8  gm/cm3,  a 
number  of  illustrative  examples  are  given  in  Table  A.l. 


4 
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TABLE  A.l.  THE  REFRACTION  COEFFICIENT  A FOR  TYPICAL  SOURCE 


a 

s 

MATERIALS 

^s 

e 

s 

A 

“bA/“o 

2.40 

1.30 

1.86 

11.53 

0.4199 

1.05 

2.80 

1.40 

1.95 

13.29 

0.4906 

1.05 

3.30 

1.90 

2.25 

15.66 

0.6130 

1.11 

4.00 

2.50 

2.40 

19.26 

0.7273 

1.09 

4.50 

2.65 

2.60 

21.71 

0.8211 

1.09 

5.10 

2.75 

2.70 

24.20 

0.9009 

1.06 

5.50 

3.10 

2.75 

26.74 

0.9475 

1.03 

For  short  period  teleseismic  body  waves  V (w)  may 
be  replaced  by  the  constant  V®.  Then  combining  all  these 
approximations  with  (A. 6) , we  find  that  teleseismic  amplitude 
is  proportional  to  c i^®.  That  is,  for  the  two  layered  earth 
model  we  deduce  that 


It  should  be  emphasized  that  it  is  the  direct  P wave 
that  scales  with  a^®.  Extending  this  to  scale  m^  implicitly 
assumes  that  the  influence  of  the  free  surface  reflected  phase 
and  details  of  the  near  source  geology  are  invarient  between 
events.  The  usefulness  of  (A. 9)  is  degraded  to  the  extent 
that  this  is  not  true. 
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Calculation  of  the  RDP 

The  RDP  is  calculated  in  the  following  manner.  If  we 
assume  that  the  displacement  is  linear  within  a given  time 
intu'*val  defined  by 


t < t < t 
1 2 


then  Eq.  (A. 2)  gives 


*(t2) 


•Ri[v\  e*]  + R*(b  - f*)  [l-«X]  + ex»(Ti)  (A. 10) 


where 


U(T1#R)  - U(TJ#R) 

a T - T 

1 2 


-U(T  , R)  X + u(X  ,R)X 


X - X 
I 2 


a {x  - x ) 

Si  j 


At  late  times  when  the  displacement  becomes  constant, 
x can  be  much  larger  than  x so  that 

2 i 


e*  * 0 


a -*■  0 


b ■ u(®,R) 


The  steady-state  value  of  the  RDP  becomes 
4<(oo)  . r2  u(»,R) 


(A. 11) 


W 
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This  result  also  follows  directly  from  Eq.  (A. 2)  since 

^(o)  = 0 (A.] 

In  an  entirely  analogous  manner,  the  reduced  velocity 
potential  (RVP)  may  be  computed  from  the  velocity-time  history 
of  a particle.  If  u(t,R)  is  radial  particle  velocity,  then 


u (t,R) 


3 y(x)  . V(t)  + V(t) 
3R  R r2  asR 


where  4'(t)  is  the  RVP. 

If  H'(t)  replaces  'Mi)  and  u(t,R)  replaces  u(t,R) 
in  Eq.  (A. 10),  then  the  algorithm  for  calculating  the  RVP 
follows  directly. 

It  is  highly  desirable  to  transform  the  equivalent 
source  into  the  frequency  domain  in  order  to  identify  the 

frequency  dependence  of  the  spectrum  required  for  equivalent 

• • 

source  comparison.  Since  ¥(»)  * 0,  then  ¥ is  bounded  at 
zero  frequency  and  the  Fourier  transform  of  $(t)  may  be  done 
numerically. 

An  important  relation  is  that  the  zero  frequency  limit 
of  the  RVP  is  equal  to  the  steady  state  value  of  the  RDP . 

Since 


w 

$ (f=0)  = f $ 


(t)  di 


f dT(i] 

J di 


Therefore,  since  ¥(0)  = 0, 


<P(0)  « VM 


Reduced  velocity  potentials  have  been  computed  for  a 
variety  of  near  source  geologic  environments.  These  calcula- 
tions were  performed  on  a Lagrangian  finite  difference  code. 
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SKIPPER  (see  Appendix  B) , that  simulates  a propagating  stress 
wave  in  one  space  dimension.  The  difference  equations  in  the 
code  are  identical  to  those  given  by  Cherry  and  Petersen. ^ 

The  code  is  capable  of  carrying  the  propagating  stress  wave 
into  the  small  displacement,  elastic  regime  and  yet  flexible 
enough  to  permit  appropriate  material  response  formulations 
in  the  large  displacement,  nonlinear  regime. 

A 

Figure  A. 2 shows  l^l  versus  frequency  for  a nuclear 
explosion  of  20  kt.  Each  source  function  in  the  figure  is 
identified  by  an  integer  and  corresponds  to  a 20  kt  explosion 
detonated  in  a geologic  environment  with  a specific  set  of 
material  properties  given  in  Section  III.  For  all  practical 
purposes  the  spectra  are  flat  within  the  teleseismic 
frequency  band  (f  < 3 Hz)  for  yields  up  to  at  least  200  kt. 

This  implies  that  teleseismic  body  wave  and  surface  wave  ampli- 
tudes should  be  directly  proportional  to  the  yield  of  the  ex- 
plosion. Deviation  from  this  linear  scaling  should  occur 
within  the  yield  interval  200  kt  < W < 1,000  kt,  depending  on 
the  source  function  used  to  represent  the  rock  material  near 
the  source. 


RVP  Amplitude  (m 
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APPENDIX  B 
THE  SKIPPER  CODE 


B.l  DESCRIPTION  OF  THE  CODE 

SKIPPER  is  a Lagrangian,  finite  difference  code  that 
simulates  a propagating  stress  wave  in  one  space  dimension. 
The  code  is  capable  of  carrying  the  propagating  stress  wave 
into  the  small  displacement,  elastic  regime  and  yet  flexible 
enough  to  permit  appropriate  material  response  formulations 
in  the  large  displacement,  nonlinear  regime. 

A schematic  of  one  computational  cycle  is  shown  in 
Fig.  B.l.  The  equation  of  motion  for  spherical  symmetry  is 


(B.l) 


where  p is  density,  u is  particle  velocity,  u is  particle 

acceleration  and  0 and  0 are  the  radial  and  tangential 

1 2 

components  of  the  stress  tensor  with  compression  positive. 

Strain  rates  in  the  radial  and  tangential  direction  are 


£ 

1 


(B.  2) 


e 

2 


(B.3 ) 


The  stress  tensor  may  be  decomposed  into  an  isotropic  com- 
ponent (P)  and  a deviator ic  component  (S^  and  S^)  where 

°i  * p - si 


i 

t 


Due  to  symmetry 

S - -2S 
1 2 


(B.  5) 
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Figure  B.l.  Cycle  of  interactions  treated  in  cal- 
culating stress  wave  propagation.  Taken 
from  Ref.  5. 
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Hooke's  law  gives 


v 

7 


S - y u (e  - e ) 

1 3 1 2 

where  k and  y are  the  bulk  and  shea:  moduli,  the  dot 
indicates  a time  derivative  along  a particle  path  and  from 
conservation  of  mass  the  rate  of  change  of  volume  is 

XT  • • 

- - e + 2e 
v 1 2 

i 

For  large  and/or  energy  dependent  deformations  Eq. 
(B. 6)  is  replaced  by 

P - f ^ (v,e)  or  P - g(n,e) 

where  v is  specific  volume. 

The  specific  internal  energy  (e)  is  obtained  from 


pe  ■ -P  - + S e + 2S  c 
v xx  12 


and  the  total  volumetric  strain  (n)  is  defined  to  be 


— v°  - v 

n ■ 

v# 


1 - — - 1 ! 

V P 

0 


The  difference  equation  used  in  the  code  that  corresponds  to 
Eqs.  (B. 1) , (B.2) , (B.3),  (B.10)  and  (B. 11)  are  the  same  as 

[7] 

those  used  by  Cherry  and  Petersen. 

Basically,  Eqs.  (B.7)  and  (B.9)  should  be  regarded  as 
an  initial  attempt  by  the  code  to  determine  P and  at 

the  current  time  cycle.  Modification  of  these  stresses  then 
depend  on  tie  constitutive  relation  used  in  the  equation  of 
state  section  of  the  code. 
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An  artificial  viscosity  is  used  in  the  momentum  and 
energy  equations  [(B.l)  and  (B.10))  and  has  the  same  form 

[Q] 

as  that  reported  by  Wilkins. 1 J This  quantity  replaces 
discontinuous  jumps  in  particle  velocity  in  the  calculation 
with  a rapid  but  continuous  change. 

B.2  SKIPPER  COMPARISON  TO  AN  ANALYTIC  SOLUTION 


Normally,  the  only  way  to  evaluate  a numerical  procedure 
is  to  examine  its  ability  to  produce  the  answer  to  a previously 
solved  problem.  Also,  in  order  to  calculate  the  equivalent 
elastic  source  function  from  a contained  nuclear  explosion, 
the  calculations  must  be  continued  well  into  the  small  dis- 
placement, elastic  regime.  In  order  for  the  calculation  to 
be  meaningful  in  this  regime  the  numerics  in  the  code  must 
not  introduce  errors  in  the  solution  when  displacements  are 
small. 

A well-known  analytic  solution,  commonly  called  Blake's 
rgi 

solution,  to  an  elastic  wave  propagation  problem  has  been 
used  for  our  analytic  comparison.  The  followina  parameters 
were  used  in  the  calculation: 


Cavity  Radius 
Density 

Compressional  Velocity 


Bulk  Modulus 


Shea’.  Modulus 


Cavity  Load  (kbar) 
Zone  Size 


ri  10  meters 


2 g/cc 
5 m/msec 


3J3.33  kbar 


125  kbar 
e (t  in  msec) 


10  cm 


Figures  B.2  and  B.3  show  particle  velocity  versus  radius  at 
1.32  msec  and  6.81  msec.  Figures  B.4  and  B.5  show  displace- 
ment and  mean  stress  versus  time  10  m from  the  cavity  wall. 
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The  numerical  results  from  SKIPPER  agree  quite  well  with  the 
analytic  solution.  Displacement  time  histories  of  the  type 
shown  in  Fig.  B.4  provide  the  input  for  the  RDP  calculation. 


r 
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